{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# Power Iteration and its Variants"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 104,
      "metadata": {},
      "outputs": [],
      "source": [
        "import numpy as np\n",
        "import numpy.linalg as la\n",
        "np.set_printoptions(precision=3)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's  prepare a matrix with some random or deliberately chosen eigenvalues:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 105,
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "[-2.668 -0.958 -0.33  -0.292 -0.186 -0.144]\n"
          ]
        }
      ],
      "source": [
        "n = 6\n",
        "\n",
        "if 1:\n",
        "    np.random.seed(70)\n",
        "    eigvecs = np.random.randn(n, n)\n",
        "    eigvals = np.sort(np.random.randn(n))\n",
        "    # Uncomment for near-duplicate largest-magnitude eigenvalue\n",
        "    # eigvals[1] = eigvals[0] + 1e-3\n",
        "\n",
        "    A = eigvecs.dot(np.diag(eigvals)).dot(la.inv(eigvecs))\n",
        "    print(eigvals)\n",
        "    \n",
        "else:\n",
        "    # Complex eigenvalues\n",
        "    np.random.seed(40)\n",
        "    A = np.random.randn(n, n)\n",
        "    print(la.eig(A)[0])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's also pick an initial vector:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 106,
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([ 2.269,  0.664,  0.899, -0.366,  0.463,  0.08 ])"
            ]
          },
          "execution_count": 106,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "x0 = np.random.randn(n)\n",
        "x0"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Power iteration"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 107,
      "metadata": {},
      "outputs": [],
      "source": [
        "x = x0"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now implement plain power iteration."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 127,
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "[ 1.329e+12 -4.740e+12 -7.728e+11  1.018e+12 -8.466e+11 -2.213e+12]\n",
            "[-3.545e+12  1.264e+13  2.062e+12 -2.716e+12  2.258e+12  5.903e+12]\n",
            "[ 9.457e+12 -3.373e+13 -5.500e+12  7.245e+12 -6.025e+12 -1.575e+13]\n",
            "[-2.523e+13  8.998e+13  1.467e+13 -1.933e+13  1.607e+13  4.201e+13]\n",
            "[ 6.730e+13 -2.400e+14 -3.914e+13  5.156e+13 -4.287e+13 -1.121e+14]\n",
            "[-1.795e+14  6.403e+14  1.044e+14 -1.375e+14  1.144e+14  2.989e+14]\n",
            "[ 4.789e+14 -1.708e+15 -2.785e+14  3.669e+14 -3.051e+14 -7.975e+14]\n",
            "[-1.278e+15  4.557e+15  7.430e+14 -9.788e+14  8.139e+14  2.127e+15]\n",
            "[ 3.408e+15 -1.216e+16 -1.982e+15  2.611e+15 -2.171e+15 -5.675e+15]\n",
            "[-9.092e+15  3.243e+16  5.288e+15 -6.965e+15  5.792e+15  1.514e+16]\n",
            "[ 2.425e+16 -8.650e+16 -1.411e+16  1.858e+16 -1.545e+16 -4.039e+16]\n",
            "[-6.470e+16  2.308e+17  3.763e+16 -4.957e+16  4.122e+16  1.077e+17]\n",
            "[ 1.726e+17 -6.156e+17 -1.004e+17  1.322e+17 -1.100e+17 -2.874e+17]\n",
            "[-4.604e+17  1.642e+18  2.678e+17 -3.527e+17  2.933e+17  7.667e+17]\n",
            "[ 1.228e+18 -4.381e+18 -7.143e+17  9.410e+17 -7.825e+17 -2.045e+18]\n",
            "[-3.277e+18  1.169e+19  1.906e+18 -2.510e+18  2.087e+18  5.456e+18]\n",
            "[ 8.741e+18 -3.117e+19 -5.083e+18  6.696e+18 -5.569e+18 -1.455e+19]\n",
            "[-2.332e+19  8.316e+19  1.356e+19 -1.786e+19  1.486e+19  3.883e+19]\n",
            "[ 6.220e+19 -2.219e+20 -3.618e+19  4.765e+19 -3.963e+19 -1.036e+20]\n",
            "[-1.659e+20  5.918e+20  9.650e+19 -1.271e+20  1.057e+20  2.763e+20]\n"
          ]
        }
      ],
      "source": [
        "for i in range(20):\n",
        "    x = A @ x\n",
        "    print(x)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "* What's the problem with this method?\n",
        "* Does anything useful come of this?\n",
        "* How do we fix it?"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Normalized power iteration"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Back to the beginning: Reset to the initial vector."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 128,
      "metadata": {},
      "outputs": [],
      "source": [
        "x = x0/la.norm(x0)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Implement normalized power iteration."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 129,
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "[-0.285  0.819  0.106 -0.172  0.149  0.431]\n",
            "[ 0.243 -0.842 -0.13   0.173 -0.156 -0.402]\n",
            "[-0.238  0.845  0.135 -0.177  0.153  0.396]\n",
            "[ 0.237 -0.845 -0.137  0.18  -0.152 -0.395]\n",
            "[-0.237  0.845  0.137 -0.181  0.151  0.395]\n",
            "[ 0.237 -0.845 -0.138  0.181 -0.151 -0.394]\n",
            "[-0.237  0.845  0.138 -0.181  0.151  0.394]\n",
            "[ 0.237 -0.845 -0.138  0.181 -0.151 -0.394]\n",
            "[-0.237  0.845  0.138 -0.181  0.151  0.394]\n",
            "[ 0.237 -0.845 -0.138  0.181 -0.151 -0.394]\n",
            "2.667662268743778\n"
          ]
        }
      ],
      "source": [
        "for i in range(10):\n",
        "    x = A @ x\n",
        "    nrm = la.norm(x)\n",
        "    x = x/nrm\n",
        "    print(x)\n",
        "\n",
        "print(nrm)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "* What do you observe about the norm?\n",
        "* What about the sign?\n",
        "* What is the vector $x$ now?\n",
        "\n",
        "Extensions:\n",
        "\n",
        "* Now try the matrix variants above.\n",
        "* Suggest a better way of estimating the eigenvalue. [Hint](https://en.wikipedia.org/wiki/Rayleigh_quotient)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "------\n",
        "### Inverse iteration\n",
        "\n",
        "What if we want the smallest eigenvalue (by magnitude)?\n",
        "\n",
        "Once again, reset to the beginning."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 143,
      "metadata": {},
      "outputs": [],
      "source": [
        "x = x0/la.norm(x0)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 144,
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "[ 0.344 -0.55  -0.094  0.23  -0.05  -0.718]\n",
            "[-0.527  0.236  0.046 -0.123 -0.066  0.803]\n",
            "[ 0.599 -0.01  -0.037 -0.031  0.127 -0.789]\n",
            "[-0.617 -0.128  0.043  0.173 -0.153  0.74 ]\n",
            "[ 0.611  0.213 -0.052 -0.284  0.161 -0.687]\n",
            "[-0.597 -0.265  0.061  0.365 -0.161  0.64 ]\n",
            "[ 0.583  0.299 -0.067 -0.424  0.16  -0.601]\n",
            "[-0.57  -0.322  0.072  0.465 -0.157  0.571]\n",
            "[ 0.558  0.337 -0.076 -0.496  0.155 -0.547]\n",
            "[-0.55  -0.348  0.078  0.517 -0.153  0.529]\n",
            "[ 0.543  0.356 -0.08  -0.533  0.151 -0.515]\n",
            "[-0.537 -0.362  0.082  0.545 -0.15   0.504]\n",
            "[ 0.533  0.366 -0.083 -0.554  0.149 -0.496]\n",
            "[-0.529 -0.369  0.084  0.561 -0.148  0.49 ]\n",
            "[ 0.527  0.371 -0.084 -0.566  0.147 -0.485]\n",
            "[-0.525 -0.373  0.085  0.57  -0.147  0.481]\n",
            "[ 0.523  0.375 -0.085 -0.573  0.146 -0.478]\n",
            "[-0.522 -0.376  0.086  0.575 -0.146  0.476]\n",
            "[ 0.521  0.376 -0.086 -0.577  0.146 -0.474]\n",
            "[-0.521 -0.377  0.086  0.578 -0.146  0.473]\n",
            "[ 0.52   0.377 -0.086 -0.579  0.146 -0.472]\n",
            "[-0.52  -0.378  0.086  0.58  -0.145  0.471]\n",
            "[ 0.519  0.378 -0.086 -0.581  0.145 -0.471]\n",
            "[-0.519 -0.378  0.086  0.581 -0.145  0.47 ]\n",
            "[ 0.519  0.379 -0.086 -0.582  0.145 -0.47 ]\n",
            "[-0.519 -0.379  0.086  0.582 -0.145  0.47 ]\n",
            "[ 0.519  0.379 -0.086 -0.582  0.145 -0.469]\n",
            "[-0.518 -0.379  0.086  0.582 -0.145  0.469]\n",
            "[ 0.518  0.379 -0.086 -0.582  0.145 -0.469]\n",
            "[-0.518 -0.379  0.086  0.583 -0.145  0.469]\n"
          ]
        }
      ],
      "source": [
        "for i in range(30):\n",
        "    x = la.solve(A, x)\n",
        "    nrm = la.norm(x)\n",
        "    x = x/nrm\n",
        "    print(x)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "* What's the computational cost per iteration?\n",
        "* Can we make this method search for a specific eigenvalue?"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "------\n",
        "### Inverse Shift iteration\n",
        "\n",
        "What if we want the eigenvalue closest to a give value $\\sigma$?\n",
        "\n",
        "Once again, reset to the beginning."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 148,
      "metadata": {},
      "outputs": [],
      "source": [
        "x = x0/la.norm(x0)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 149,
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "[ 0.011 -0.809 -0.191  0.223 -0.186 -0.474]\n",
            "[-0.127  0.775  0.149 -0.232  0.16   0.531]\n",
            "[ 0.192 -0.728 -0.123  0.246 -0.135 -0.582]\n",
            "[-0.246  0.674  0.103 -0.261  0.108  0.628]\n",
            "[ 0.295 -0.616 -0.085  0.271 -0.08  -0.669]\n",
            "[-0.338  0.556  0.07  -0.277  0.052  0.702]\n",
            "[ 0.376 -0.497 -0.057  0.277 -0.026 -0.728]\n",
            "[-0.41   0.441  0.047 -0.273  0.002  0.749]\n",
            "[ 0.44  -0.387 -0.038  0.264  0.019 -0.765]\n",
            "[-0.465  0.338  0.031 -0.252 -0.039  0.777]\n",
            "[ 0.488 -0.291 -0.026  0.237  0.056 -0.786]\n",
            "[-0.508  0.249  0.022 -0.219 -0.071  0.792]\n",
            "[ 0.525 -0.209 -0.019  0.199  0.084 -0.796]\n",
            "[-0.54   0.172  0.017 -0.178 -0.095  0.798]\n",
            "[ 0.553 -0.138 -0.016  0.155  0.105 -0.799]\n",
            "[-0.565  0.107  0.016 -0.132 -0.114  0.799]\n",
            "[ 0.575 -0.078 -0.016  0.108  0.122 -0.798]\n",
            "[-0.584  0.05   0.016 -0.084 -0.129  0.796]\n",
            "[ 0.591 -0.025 -0.017  0.06   0.135 -0.792]\n",
            "[-0.597  0.001  0.018 -0.035 -0.14   0.789]\n",
            "[ 0.603  0.021 -0.02   0.011  0.144 -0.784]\n",
            "[-0.607 -0.041  0.022  0.013 -0.148  0.779]\n",
            "[ 0.611  0.061 -0.024 -0.036  0.151 -0.774]\n",
            "[-0.614 -0.079  0.025  0.059 -0.154  0.768]\n",
            "[ 0.616  0.096 -0.027 -0.082  0.156 -0.761]\n",
            "[-0.618 -0.112  0.03   0.103 -0.158  0.755]\n",
            "[ 0.619  0.126 -0.032 -0.124  0.16  -0.748]\n",
            "[-0.619 -0.14   0.034  0.145 -0.161  0.741]\n",
            "[ 0.619  0.154 -0.036 -0.164  0.162 -0.734]\n",
            "[-0.619 -0.166  0.038  0.183 -0.163  0.726]\n"
          ]
        }
      ],
      "source": [
        "sigma = 1\n",
        "A_sigma = A-sigma*np.eye(A.shape[0])\n",
        "for i in range(30):\n",
        "    x = la.solve(A_sigma, x)\n",
        "    nrm = la.norm(x)\n",
        "    x = x/nrm\n",
        "    print(x)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "--------------\n",
        "### Rayleigh quotient iteration\n",
        "\n",
        "Can we feed an estimate of the current approximate eigenvalue back into the calculation? (Hint: Rayleigh quotient)\n",
        "\n",
        "Reset once more."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 153,
      "metadata": {},
      "outputs": [],
      "source": [
        "x = x0/la.norm(x0)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Run this cell in-place (Ctrl-Enter) many times."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 154,
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "[ 0.063 -0.792 -0.173  0.23  -0.173 -0.505]\n",
            "[-0.191  0.726  0.126 -0.248  0.133  0.585]\n",
            "[ 0.521 -0.254 -0.051  0.121  0.058 -0.802]\n",
            "[-0.544 -0.35   0.081  0.533 -0.15   0.519]\n",
            "[ 0.52   0.378 -0.086 -0.58   0.146 -0.472]\n",
            "[-0.518 -0.379  0.086  0.583 -0.145  0.469]\n",
            "[ 0.518  0.379 -0.086 -0.583  0.145 -0.469]\n",
            "[-0.518 -0.379  0.086  0.583 -0.145  0.469]\n",
            "[-0.518 -0.379  0.086  0.583 -0.145  0.469]\n",
            "[ 0.518  0.379 -0.086 -0.583  0.145 -0.469]\n"
          ]
        }
      ],
      "source": [
        "for i in range(10):\n",
        "    sigma = np.dot(x, np.dot(A, x))/np.dot(x, x)\n",
        "    x = la.solve(A-sigma*np.eye(n), x)\n",
        "    x = x/la.norm(x)\n",
        "    print(x)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "* What's a reasonable stopping criterion?\n",
        "* Computational downside of this iteration?"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.6.5"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 1
}